sep_smooths_year2
Load all required libraries.
library(tidyverse)## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --## v ggplot2 3.3.5 v purrr 0.3.4
## v tibble 3.1.3 v dplyr 1.0.7
## v tidyr 1.1.3 v stringr 1.4.0
## v readr 2.0.0 v forcats 0.5.1## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()library(plotly)##
## Attaching package: 'plotly'## The following object is masked from 'package:ggplot2':
##
## last_plot## The following object is masked from 'package:stats':
##
## filter## The following object is masked from 'package:graphics':
##
## layoutlibrary(broom)Read in raw data from RDS.
raw_data <- readRDS("./year2.RDS")Make a few small modifications to names and data for visualizations.
final_data <- raw_data %>% mutate(log_copy_per_L = log10(mean_copy_num_L)) %>%
rename(Facility = wrf) %>%
mutate(Facility = recode(Facility,
"NO" = "WRF A",
"MI" = "WRF B",
"CC" = "WRF C"))Seperate the data by gene target to ease layering in the final plot
#make three data layers
only_positives <<- subset(final_data, (!is.na(final_data$Facility)))
only_n1 <- subset(only_positives, target == "N1")
only_n2 <- subset(only_positives, target == "N2")
only_background <<-final_data %>%
select(c(date, cases_cum_clarke, new_cases_clarke, X7_day_ave_clarke)) %>%
group_by(date) %>% summarise_if(is.numeric, mean)
#specify fun colors
background_color <- "#7570B3"
seven_day_ave_color <- "#E6AB02"
marker_colors <- c("N1" = '#1B9E77',"N2" ='#D95F02')#remove facilty C for now
#only_n1 <- only_n1[!(only_n1$Facility == "WRF C"),]
#only_n2 <- only_n2[!(only_n2$Facility == "WRF C"),]
only_n1 <- only_n1[!(only_n1$Facility == "WRF A" & only_n1$date == "2020-11-02"), ]
only_n2 <- only_n2[!(only_n2$Facility == "WRF A" & only_n2$date == "2020-11-02"), ]Build the main plot
#first layer is the background epidemic curve
p1 <- only_background %>%
plotly::plot_ly() %>%
plotly::add_trace(x = ~date, y = ~new_cases_clarke,
type = "bar",
hoverinfo = "text",
text = ~paste('</br> Date: ', date,
'</br> Daily Cases: ', new_cases_clarke),
alpha = 0.5,
name = "Daily Reported Cases",
color = background_color,
colors = background_color,
showlegend = FALSE) %>%
layout(yaxis = list(title = "Clarke County Daily Cases", showline=TRUE)) %>%
layout(legend = list(orientation = "h", x = 0.2, y = -0.3))
#renders the main plot layer two as seven day moving average
p1 <- p1 %>% plotly::add_trace(x = ~date, y = ~X7_day_ave_clarke,
type = "scatter",
mode = "lines",
hoverinfo = "text",
text = ~paste('</br> Date: ', date,
'</br> Seven-Day Moving Average: ', X7_day_ave_clarke),
name = "Seven Day Moving Average Athens",
line = list(color = seven_day_ave_color),
showlegend = FALSE)
#renders the main plot layer three as positive target hits
p2 <- plotly::plot_ly() %>%
plotly::add_trace(x = ~date, y = ~mean_copy_num_L,
type = "scatter",
mode = "markers",
hoverinfo = "text",
text = ~paste('</br> Date: ', date,
'</br> Facility: ', Facility,
'</br> Target: ', target,
'</br> Copies/L: ', round(mean_copy_num_L, digits = 2)),
data = only_n1,
symbol = ~Facility,
marker = list(color = '#1B9E77', size = 8, opacity = 0.65),
showlegend = FALSE) %>%
plotly::add_trace(x = ~date, y = ~mean_copy_num_L,
type = "scatter",
mode = "markers",
hoverinfo = "text",
text = ~paste('</br> Date: ', date,
'</br> Facility: ', Facility,
'</br> Target: ', target,
'</br> Copies/L: ', round(mean_copy_num_L, digits = 2)),
data = only_n2,
symbol = ~Facility,
marker = list(color = '#D95F02', size = 8, opacity = 0.65),
showlegend = FALSE) %>%
layout(yaxis = list(title = "SARS CoV-2 Copies/L",
showline = TRUE,
type = "log",
dtick = 1,
automargin = TRUE)) %>%
layout(legend = list(orientation = "h", x = 0.2, y = -0.3))
#adds the limit of detection dashed line
p2 <- p2 %>% plotly::add_segments(x = as.Date("2021-06-30"),
xend = ~max(date + 10),
y = 3571.429, yend = 3571.429,
opacity = 0.35,
line = list(color = "black", dash = "dash")) %>%
layout(annotations = list(x = as.Date("2021-06-30"), y = 3.8, xref = "x", yref = "y",
text = "Limit of Detection", showarrow = FALSE))
p1 p2Combine the two main plot pieces as a subplot
#seperate n1 and n2 frames by site
#n1
wrf_a_only_n1 <- subset(only_n1, Facility == "WRF A")
wrf_b_only_n1 <- subset(only_n1, Facility == "WRF B")
wrf_c_only_n1 <- subset(only_n1, Facility == "WRF C")
#n2
wrf_a_only_n2 <- subset(only_n2, Facility == "WRF A")
wrf_b_only_n2 <- subset(only_n2, Facility == "WRF B")
wrf_c_only_n2 <- subset(only_n2, Facility == "WRF C")
#rejoin the old data frames then seperate in to averages for each plant.
wrfa_both <- full_join(wrf_a_only_n1, wrf_a_only_n2)%>%
select(c(date, mean_total_copies)) %>%
group_by(date) %>%
summarize_if(is.numeric, mean) %>%
ungroup() %>%
mutate(log_total_copies_both = log10(mean_total_copies))## Joining, by = c("date", "new_cases_clarke", "cases_cum_clarke", "X7_day_ave_clarke", "Facility", "collection_num", "target", "mean_copy_num_uL_rxn", "mean_copy_num_L", "sd_L", "mean_total_copies", "sd_total_copies", "log_copy_per_L")wrfb_both <- full_join(wrf_b_only_n1, wrf_b_only_n2)%>%
select(c(date, mean_total_copies)) %>%
group_by(date) %>%
summarize_if(is.numeric, mean) %>%
ungroup() %>%
mutate(log_total_copies_both = log10(mean_total_copies))## Joining, by = c("date", "new_cases_clarke", "cases_cum_clarke", "X7_day_ave_clarke", "Facility", "collection_num", "target", "mean_copy_num_uL_rxn", "mean_copy_num_L", "sd_L", "mean_total_copies", "sd_total_copies", "log_copy_per_L")wrfc_both <- full_join(wrf_c_only_n1, wrf_c_only_n2)%>%
select(c(date, mean_total_copies)) %>%
group_by(date) %>%
summarize_if(is.numeric, mean) %>%
ungroup() %>%
mutate(log_total_copies_both = log10(mean_total_copies))## Joining, by = c("date", "new_cases_clarke", "cases_cum_clarke", "X7_day_ave_clarke", "Facility", "collection_num", "target", "mean_copy_num_uL_rxn", "mean_copy_num_L", "sd_L", "mean_total_copies", "sd_total_copies", "log_copy_per_L")#get max date
maxdate <- max(wrfa_both$date)
mindate <- min(wrfa_both$date)Build loess smoothing figures figures
This makes the individual plots
#**************************************WRF A PLOT**********************************************
#add trendlines
#extract data from geom_smooth
#both extract
# *********************************span 0.6***********************************
#*****************Must always update the n = TOTAL NUMBER OF DAYS*************************
extract_botha <- ggplot(wrfa_both, aes(x = date, y = log_total_copies_both)) +
stat_smooth(aes(outfit=fit_botha<<-..y..), method = "loess", color = '#1B9E77',
span = 0.3, n = 99)## Warning: Ignoring unknown aesthetics: outfit#look at the fits to align dates and total observations
#both
extract_botha## `geom_smooth()` using formula 'y ~ x'
fit_botha## [1] 11.58564 11.60000 11.61801 11.63926 11.66330 11.68971 11.71806 11.75010
## [9] 11.78702 11.82757 11.87049 11.91453 11.95843 12.03249 12.10724 12.13889
## [17] 12.15115 12.15428 12.15857 12.17429 12.21455 12.25773 12.28978 12.32368
## [25] 12.35808 12.39165 12.42303 12.44251 12.46303 12.50194 12.55031 12.60132
## [33] 12.64814 12.68397 12.71481 12.74188 12.75600 12.76090 12.76200 12.76476
## [41] 12.77462 12.77183 12.77153 12.79412 12.82114 12.85209 12.88648 12.92382
## [49] 12.99078 13.05858 13.10552 13.16032 13.21492 13.26123 13.29118 13.30298
## [57] 13.30677 13.29799 13.27154 13.23642 13.20158 13.17600 13.15873 13.14172
## [65] 13.12575 13.11164 13.09790 13.08303 13.06553 13.04390 13.01052 12.98972
## [73] 12.97220 12.95495 12.93493 12.90910 12.88436 12.85639 12.81664 12.77136
## [81] 12.72291 12.67364 12.62591 12.56363 12.50171 12.45225 12.39850 12.34604
## [89] 12.30040 12.26713 12.25462 12.24651 12.23295 12.22430 12.21934 12.21687
## [97] 12.21567 12.21542 12.21710#assign fits to a vector
both_trenda <- fit_botha
#extract y min and max for each
limits_botha <- ggplot_build(extract_botha)$data## `geom_smooth()` using formula 'y ~ x'limits_botha <- as.data.frame(limits_botha)
both_ymina <- limits_botha$ymin
both_ymaxa <- limits_botha$ymax
#reassign dataframes (just to be safe)
work_botha <- wrfa_both
#fill in missing dates to smooth fits
work_botha <- work_botha %>% complete(date = seq(min(date), max(date), by = "1 day"))
date_vec_botha <- work_botha$date
#create a new smooth dataframe to layer
smooth_frame_botha <- data.frame(date_vec_botha, both_trenda, both_ymina, both_ymaxa)#WRF A
#plot smooth frames
p_wrf_a <- plotly::plot_ly() %>%
plotly::add_lines(x = ~date_vec_botha, y = ~both_trenda,
data = smooth_frame_botha,
hoverinfo = "text",
text = ~paste('</br> Date: ', date_vec_botha,
'</br> Median Log Copies: ', round(both_trenda, digits = 2)),
line = list(color = '#1B9E77', size = 8, opacity = 0.65),
showlegend = FALSE) %>%
layout(xaxis = list(range = c(mindate - 7, maxdate + 7))) %>% #buffer here
plotly::add_ribbons(x ~date_vec_botha, ymin = ~both_ymina, ymax = ~both_ymaxa,
showlegend = FALSE,
opacity = 0.25,
hoverinfo = "text",
text = ~paste('</br> Date: ', date_vec_botha, #leaving in case we want to change
'</br> Max Log Copies: ', round(both_ymaxa, digits = 2),
'</br> Min Log Copies: ', round(both_ymina, digits = 2)),
name = "",
fillcolor = '#1B9E77',
line = list(color = '#1B9E77')) %>%
layout(yaxis = list(title = "Total Log10 SARS CoV-2 Copies",
showline = TRUE,
automargin = TRUE)) %>%
layout(xaxis = list(title = "Date")) %>%
layout(title = "WRF A") %>%
plotly::add_markers(x = ~date, y = ~log_total_copies_both,
data = wrfa_both,
hoverinfo = "text",
showlegend = FALSE,
text = ~paste('</br> Date: ', date,
'</br> Actual Log Copies: ', round(log_total_copies_both, digits = 2)),
marker = list(color = '#1B9E77', size = 6, opacity = 0.65))
p_wrf_asave(p_wrf_a, file = "./site_objects/wrf_a_year2.rda")#**************************************WRF B PLOT**********************************************
#add trendlines
#extract data from geom_smooth
#both extract
# *********************************span 0.6***********************************
#*****************Must always update the n = TOTAL NUMBER OF DAYS*************************
extract_bothb <- ggplot(wrfb_both, aes(x = date, y = log_total_copies_both)) +
stat_smooth(aes(outfit=fit_bothb<<-..y..), method = "loess", color = '#D95F02',
span = 0.3, n = 99)## Warning: Ignoring unknown aesthetics: outfit#look at the fits to align dates and total observations
#both
extract_bothb## `geom_smooth()` using formula 'y ~ x'
fit_bothb## [1] 10.65612 10.78161 10.89978 11.01104 11.11579 11.21444 11.30738 11.39240
## [9] 11.46823 11.53666 11.59950 11.65856 11.71564 11.74419 11.76780 11.81324
## [17] 11.85758 11.90334 11.95303 12.00917 12.07917 12.15363 12.23462 12.33189
## [25] 12.43421 12.53037 12.60914 12.68551 12.75422 12.79818 12.83009 12.85311
## [33] 12.87040 12.88513 12.88333 12.87701 12.87381 12.85816 12.83741 12.81891
## [41] 12.81001 12.80345 12.80008 12.81057 12.83007 12.85485 12.88122 12.90544
## [49] 12.93803 12.96965 12.99129 13.01441 13.03626 13.05405 13.06504 13.07465
## [57] 13.08024 13.07114 13.05053 13.02444 12.99888 12.97988 12.96328 12.94714
## [65] 12.93915 12.93816 12.94018 12.94121 12.93729 12.92440 12.87877 12.85182
## [73] 12.81391 12.77469 12.74382 12.73092 12.74658 12.76840 12.79564 12.84455
## [81] 12.90078 12.94997 12.97778 12.98263 12.97995 12.97016 12.94852 12.91849
## [89] 12.88353 12.84710 12.80859 12.76566 12.71855 12.66699 12.61113 12.55113
## [97] 12.48715 12.41913 12.34683#assign fits to a vector
both_trendb <- fit_bothb
#extract y min and max for each
limits_bothb <- ggplot_build(extract_bothb)$data## `geom_smooth()` using formula 'y ~ x'limits_bothb <- as.data.frame(limits_bothb)
both_yminb <- limits_bothb$ymin
both_ymaxb <- limits_bothb$ymax
#reassign dataframes (just to be safe)
work_bothb <- wrfb_both
#fill in missing dates to smooth fits
work_bothb <- work_bothb %>% complete(date = seq(min(date), max(date), by = "1 day"))
date_vec_bothb <- work_bothb$date
#create a new smooth dataframe to layer
smooth_frame_bothb <- data.frame(date_vec_bothb, both_trendb, both_yminb, both_ymaxb)#WRF B
#plot smooth frames
p_wrf_b <- plotly::plot_ly() %>%
plotly::add_lines(x = ~date_vec_bothb, y = ~both_trendb,
data = smooth_frame_bothb,
hoverinfo = "text",
text = ~paste('</br> Date: ', date_vec_bothb,
'</br> Median Log Copies: ', round(both_trendb, digits = 2)),
line = list(color = '#D95F02', size = 8, opacity = 0.65),
showlegend = FALSE) %>%
layout(xaxis = list(range = c(mindate - 7, maxdate + 7))) %>% #buffer here
plotly::add_ribbons(x ~date_vec_bothb, ymin = ~both_yminb, ymax = ~both_ymaxb,
showlegend = FALSE,
opacity = 0.25,
hoverinfo = "text",
text = ~paste('</br> Date: ', date_vec_bothb, #leaving in case we want to change
'</br> Max Log Copies: ', round(both_ymaxb, digits = 2),
'</br> Min Log Copies: ', round(both_yminb, digits = 2)),
name = "",
fillcolor = '#D95F02',
line = list(color = '#D95F02')) %>%
layout(yaxis = list(title = "Total Log10 SARS CoV-2 Copies",
showline = TRUE,
automargin = TRUE)) %>%
layout(xaxis = list(title = "Date")) %>%
layout(title = "WRF B") %>%
plotly::add_markers(x = ~date, y = ~log_total_copies_both,
data = wrfb_both,
hoverinfo = "text",
showlegend = FALSE,
text = ~paste('</br> Date: ', date,
'</br> Actual Log Copies: ', round(log_total_copies_both, digits = 2)),
marker = list(color = '#D95F02', size = 6, opacity = 0.65))
p_wrf_bsave(p_wrf_b, file = "./site_objects/wrf_b_year2.rda")#**************************************WRF C PLOT********************************************** #add trendlines #extract data from geom_smooth # *********************************span 0.6*********************************** #*****************Must always update the n = TOTAL NUMBER OF DAYS*************************
extract_bothc <- ggplot(wrfc_both, aes(x = date, y = log_total_copies_both)) +
stat_smooth(aes(outfit=fit_bothc<<-..y..), method = "loess", color = '#E7298A',
span = 0.3, n = 99)## Warning: Ignoring unknown aesthetics: outfit#look at the fits to align dates and total observations
#both
extract_bothc## `geom_smooth()` using formula 'y ~ x'
fit_bothc## [1] 10.58009 10.71541 10.83865 10.95133 11.05495 11.15102 11.24105 11.31978
## [9] 11.38366 11.43675 11.48313 11.52687 11.57204 11.53946 11.50291 11.58318
## [17] 11.72118 11.88046 12.02458 12.11711 12.15209 12.17113 12.16537 12.11836
## [25] 12.05367 11.99487 11.96554 11.92634 11.89422 11.91511 11.95748 12.01316
## [33] 12.07396 12.13169 12.23020 12.32884 12.38862 12.45425 12.51644 12.56591
## [41] 12.59340 12.57960 12.55775 12.53208 12.47482 12.40742 12.35133 12.32799
## [49] 12.33269 12.34731 12.39294 12.47480 12.57043 12.65739 12.71324 12.75593
## [57] 12.78990 12.79020 12.76853 12.73716 12.70840 12.69451 12.69222 12.69001
## [65] 12.69288 12.70273 12.71497 12.72501 12.72826 12.72013 12.69000 12.66397
## [73] 12.62730 12.58514 12.54268 12.50509 12.45702 12.41052 12.38343 12.36184
## [81] 12.34433 12.32952 12.31598 12.30833 12.30188 12.29154 12.28255 12.27501
## [89] 12.26904 12.26472 12.25809 12.25331 12.25419 12.25663 12.26116 12.26828
## [97] 12.27853 12.29213 12.30875#assign fits to a vector
both_trendc <- fit_bothc
#extract y min and max for each
limits_bothc <- ggplot_build(extract_bothc)$data## `geom_smooth()` using formula 'y ~ x'limits_bothc <- as.data.frame(limits_bothc)
both_yminc <- limits_bothc$ymin
both_ymaxc <- limits_bothc$ymax
#reassign dataframes (just to be safe)
work_bothc <- wrfc_both
#fill in missing dates to smooth fits
work_bothc <- work_bothc %>% complete(date = seq(min(date), max(date), by = "1 day"))
date_vec_bothc <- work_bothc$date
#create a new smooth dataframe to layer
smooth_frame_bothc <- data.frame(date_vec_bothc, both_trendc, both_yminc, both_ymaxc)#WRF C
#plot smooth frames
p_wrf_c <- plotly::plot_ly() %>%
plotly::add_lines(x = ~date_vec_bothc, y = ~both_trendc,
data = smooth_frame_bothc,
hoverinfo = "text",
text = ~paste('</br> Date: ', date_vec_bothc,
'</br> Median Log Copies: ', round(both_trendc, digits = 2)),
line = list(color = '#E7298A', size = 8, opacity = 0.65),
showlegend = FALSE) %>%
layout(xaxis = list(range = c(mindate - 7, maxdate + 7))) %>% #buffer here
plotly::add_ribbons(x ~date_vec_bothc, ymin = ~both_yminc, ymax = ~both_ymaxc,
showlegend = FALSE,
opacity = 0.25,
hoverinfo = "text",
text = ~paste('</br> Date: ', date_vec_bothc, #leaving in case we want to change
'</br> Max Log Copies: ', round(both_ymaxc, digits = 2),
'</br> Min Log Copies: ', round(both_yminc, digits = 2)),
name = "",
fillcolor = '#E7298A',
line = list(color = '#E7298A')) %>%
layout(yaxis = list(title = "Total Log10 SARS CoV-2 Copies",
showline = TRUE,
automargin = TRUE)) %>%
layout(xaxis = list(title = "Date")) %>%
layout(title = "WRF C") %>%
plotly::add_markers(x = ~date, y = ~log_total_copies_both,
data = wrfc_both,
hoverinfo = "text",
showlegend = FALSE,
text = ~paste('</br> Date: ', date,
'</br> Actual Log Copies: ', round(log_total_copies_both, digits = 2)),
marker = list(color = '#E7298A', size = 6, opacity = 0.65))
p_wrf_csave(p_wrf_c, file = "./site_objects/wrf_c_year2.rda")keeping in case
#save(wrfa_both, file = "./plotly_objs/wrfa_both.rda")
#save(wrfb_both, file = "./plotly_objs/wrfb_both.rda")
#save(wrfc_both, file = "./plotly_objs/wrfc_both.rda")#save(date_vec_botha, file = "./plotly_objs/date_vec_botha.rda")
#save(date_vec_bothb, file = "./plotly_objs/date_vec_bothb.rda")
#save(date_vec_bothc, file = "./plotly_objs/date_vec_bothc.rda")#save(both_ymina, file = "./plotly_objs/both_ymina.rda")
#save(both_ymaxa, file = "./plotly_objs/both_ymaxa.rda")
#save(both_yminb, file = "./plotly_objs/both_yminb.rda")
#save(both_ymaxb, file = "./plotly_objs/both_ymaxb.rda")
#save(both_yminc, file = "./plotly_objs/both_yminc.rda")
#save(both_ymaxc, file = "./plotly_objs/both_ymaxc.rda")